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Abstract 

This work is concerned with approximating the smallest eigenvalue of a parameter- 
dependent Hermitian matrix A{fj,) for many parameter values fj, € The design 
of reliable and efficient algorithms for addressing this task is of importance in a 
variety of applications. Most notably, it plays a crucial role in estimating the error 
of reduced basis methods for parametrized partial differential equations. The cur¬ 
rent state-of-the-art approach, the so called Successive Constraint Method (SCM), 
addresses affine linear parameter dependencies by combining sampled Rayleigh quo¬ 
tients with linear programming techniques. In this work, we propose a subspace 
approach that additionally incorporates the sampled eigenvectors of A(fi) and im¬ 
plicitly exploits their smoothness properties. Like SCM, our approach results in 
rigorous lower and upper bounds for the smallest eigenvalues on D. Theoretical and 
experimental evidence is given to demonstrate that our approach represents a signif¬ 
icant improvement over SCM in the sense that the bounds are often much tighter, 
at negligible additional cost. 

Keywords, parameter-dependent eigenvalue problem, Hermitian matrix, subspace ac¬ 
celeration, Successive Constraint Method, quadratic residual bound 

1 Introduction 

Let A : D —>■ be a matrix-valued function on a compact subset D C such that 

A(/r) is Hermitian for every jj, ^ D. We aim at approximating the smallest eigenvalue, 

Amin(A(//)), ^ D, (1) 

of A{fj,) for many different values of fi G D. We consider a large-scale setting, where 
applying a standard eigensolver, such as the Lanczos method [2], is computationally 
feasible for a few values of /r but would become too expensive for many (e.g., thou¬ 
sand) parameter values. Guiding our developments, an important application of 0 
consists of estimating the coercivity constant for parametrized elliptic partial differen¬ 
tial equations (PDEs), see, e.g., [M]. In turn, these estimates can be used to construct 
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reliable a posteriori error estimates in the reduced basis method (RBM) for solving such 
PDEs. For more general PDFs, the coercivity constant needs to be replaced by the inf- 
sup constant, which - after discretization ~ corresponds to the smallest singular value 
of a general nonsymmetric matrix A{^) or, equivalently, to the smallest eigenvalue of 
the Hermitian matrix A[ijl)*A{^). Other applications that require the solution of such 
parameter-dependent eigenvalue and singular value problems include the computation 
of pseudospectra [371 Part IX], the method of particular solutions [Ij, and the eigenvalue 
analysis of waveguides [6]. The related problem of optimizing the smallest eigenvalue(s) 
of a parameter-dependent Hermitian matrix appears in a large variety of applications: 
One-parameter optimization problems play a critical role in the design of numerical 
methods |33j and robust control [2T]; multi-parameter optimization problems arise from 
semidefinite programming |9] and graph partitioning HSIE!- 

Without any further assumptions on the dependence of H(//) on the solution of Q 
is computationally intractable, especially when P is large. An assumption commonly 
found in RBM is that A(/x) admits an affine linear decomposition with respect to /U. 

Assumption 1.1 (Affine linear decomposition). Given Q G N, the Hermitian matrix 
A{fj.) admits a decomposition of the form 

A{iJ,) = 9i{p,)Ai + ■ ■ ■ + 6 q{p,)Aq, V/r G T), (2) 

for Hermitian matrices Ai, ..., Aq G and functions 6i,... ,6 q : H i—)• M. 

Assumption 1 1.1 1 holds with small Q for a number of important applications, including 
PDEs with parametrized coefficients on disjoint subdomains |3l]. Even when A(p,) does 
not satisfy this assumption, it may still be possible to approximate it very well by a 
short affine linear decomposition using, e.g., the Empirical Interpolation Method [2j. 

One of the simplest approaches to address Q is to use Gershgorin’s theorem [2] for 
estimating the smallest eigenvalue, but the accuracy of the resulting estimate is usually 
insufficient and limits the scope of applications severely. Within the context of RBM, 
a number of approaches have been developed that go beyond this simple estimate by 
making use of Assumption o For example, eigenvalue perturbation analysis can be 
used to locally approximate the smallest eigenvalues [28l |38] . The Successive Constraint 
Method (SCM; see [l2]) is currently the most commonly used approach within RBM, 
probably due to its generality and relative simplicity. Variants of SCM for computing 
smallest singular values can be found in |35|lll). while an extension of SCM to non-linear 
problems and alternative heuristic strategies have been proposed in |24j . 

If A depends analytically on n then the smallest eigenvalue inherits this property 
if Amin(A(//)) remains simple [15]. As shown in [I], the analyticity can be used to 
approximate Amin(A(//)) very well by high-order Legendre polynomials (for P = 1) 
or sparse tensor products of Legendre polynomials (for P > 1 if D is a hypercube). 
Requiring Amin(A(/u)) to stay simple on the whole of D is, however, a rather strong 
condition. In general, there are eigenvalue crossings at which Aniin(A(/i)) is Lipschitz 
continuous only; see m for a recently proposed eigenvalue optimization method that 
takes this piecewise regularity of Amin(A(^)) into account. For larger P, keeping track 
of eigenvalue crossings explicitly appears to be a rather daunting task and we therefore 
aim at a method for solving 0 that benefits only implicitly from piecewise regularity. 

The approach proposed in this paper can be summarized as follows. Given J parame¬ 
ter samples /r 2 , • • •, we consider the subspace V containing eigenvectors belonging 
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to one or several smallest eigenvalues of for i = 1,..., J. The smallest Ritz value 

of A{iJ.) with respect to V immediately yields an upper bound for Aniin(^(/^))- A lower 
bound is obtained by combining this upper bound with a perturbation argument. To 
apply such an argument requires, however, knowledge on the involved eigenvalue gap. 
We show that this gap can be estimated by adapting the linear programming approach 
used in SCM. The difference between upper and lower bounds constitutes the error esti¬ 
mate that drives the greedy strategy for selecting the next parameter sample /J-j+i- The 
whole procedure is stopped once the error estimate is uniformly small on D or, rather, 
on a surrogate of D. 

As we will see in the numerical experiments, our subspace approach accelerates con¬ 
vergence significantly compared to SCM. Subspace approaches based on additional con¬ 
ditions on the parameter dependencies have been proposed in |23[ IHT] . In the context of 
eigenvalue optimization problems, subspace acceleration has been discussed in mm- 

The rest of this paper is organized as follows. In Section we first give a brief 
overview of SCM. We then discuss its interpolation properties and point out a limitation 
on the quality of the lower bounds that can possibly be attained when solely using the 
information taken into account by SCM. In Section we present our novel subspace- 
accelerated approach for solving 0. Furthermore, we show that the new approach 
has better interpolation properties than SCM. Motivated by the fast convergence of the 
upper bounds in the novel approach, we also introduce residual-based lower bounds which 
are less reliable but sometimes converge much faster. In Sections and we present 
numerical experiments and discuss the application of our algorithm to the computation 
of coercivity and inf-sup constants. 


2 Successive Constraint Method 


In the following, we recall the Successive Constraint Method (SCM) from [12] and derive 
some new properties. The basic idea of SCM is to construct reduced-order models for 0 
that allow the efficient evaluation of of lower and upper bounds for Amin Being a 

consequence of Assumption |1.1[ the following characterization of the smallest eigenvalue 
is central to SCM: 


-^min(A(/r)) 


mm 

u^O 

min 

uec^ 

ii#0 


u* A{^)u 
u*u 




Q 

min 

uec-V 
u^o 9=1 

= min6l(/i)'^y, 
y&y 




U* AqU 
U*U 


( 3 ) 


where we have defined the vector-valued functions 9 : D ^ R : \ {0} —)■ 


as 


9{fi) :=[9i{n), •••, , R{u) 


u*Aiu 

u*u 


U*AqU 

u*u 



( 4 ) 


and set y := im(R). It follows from Q that the computation of Amin(A(/i)) is eqnivalent 
to optimizing a linear functional over y. The constraint set y is called the joint numerical 
range of Ai,..., Aq, which is generally not convex; see [8|. Thus standard optimization 
techniques cannot be used to reliably solve In SCM, the set y is approximated from 
above and from below by convex polyhedra. In turn, this allows for the use of linear 
programming (LP) techniques to yield lower and upper bounds. 
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2.1 Basic idea of SCM 


Given J parameter values Cj = {^i,... C D, let us suppose we have computed 
the corresponding eigenpairs (Ai,ui),..., (Aj,uj), that is, A* is the smallest eigenvalue 
of with eigenvector Vi G C'^. We now describe how SCM uses this information to 

approximate the set 3^ defined above. 

Clearly, 

3^UB(Cj):={i?(ui):f = l,...,J} (5) 

is a subset of y. Optimizing Q over 3 ^ub(C'j) instead of y thus yields an upper bound for 
Amin(^(Ai))- Note that this is equivalent to optimizing over the convex hull of 3 ^ub(C'j) 5 
since a solution of the LP can always be attained at a vertex of the convex polyhedron. 
To get a lower bound, we first dehne the bounding box 

B = [Amin (Al), Amax(^l)] X • • • X [A min (^Q),A^a.(AlQ)] CRQ. (6) 

By the minimax characterization of eigenvalues we have y C B, but this approximation 
is often too crude and we will instead work with 

3^lb(Cj) ■.= {yeB: > Xi,i = 1,..., J}. 

The property T C TubC^j) follows from the fact that every y = R{uy) G y satisfies 
9{lJ,i)'^y = u*A{fj,i)uy/u*Uy > m.muU*A{y,i)u/u*u = A*. The minimax characterization 
also implies that the convex polyhedron TubC^j) is tangential to y. 

With the sets defined above, we let 

XuB{y;Cj) := min 6{^ify, Alb(/^; C'j) := min 0{yfy. (7) 

y&yvji(Cj) y&yi.B{Cj) 

Since Tub(C'j) ^ T ^ T’lb(C'j), it follows that 

Alb(m; Cj) ^ Amin {A{y)) < Aub(/^; Cj) 

for every y G D. While the evaluation of Aub(/u; Cj) is trivial, the evaluation of 
Alb(/^; Cj) requires the solution of an LP; see Figure for an illustration. 


2.2 Error estimate and sampling strategy 


Assessing the quality of the bounds 0 on the entire, usually continuous parameter 
domain D is, in general, an infeasible task. A common strategy in SCM, we substitute D 
by a training set 3 C D that contains hnitely many (usually, a few thousand) parameter 
samples. We then measure the quality of the bounds by estimating the largest relative 
difference: 


max 

yes 


Aub(/^; Cj) - Alb(/^; Cj) 
|AuB(f^; Cj)\ 


( 8 ) 


If Q is not sufficiently small, SCM enlarges Cj by a parameter that attains the 
maximum in Q and recomputes the bounds Q . The resulting greedy sampling strategy 
is summarized in Algorithm 
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Figure 1: Illustration of the LP defining the lower bound Alb(m;C'j) for Q = 2 and 
J = 2. 


Algorithm 1 Successive Constraint Method 


Input: Training set H, affine linear decomposition such that A(/r) = 9i[^)Ai + 
9q{^)Aq is Hermitian for every ^ G H. Relative error tolerance escM- 
Output: Set Cj C H with corresponding eigenpairs {Xi,Vi), such that 

< "SCM for every p e =. 

compute Amin(^q), Amax(Aq) foT q = 1,..., Q, defining B according to Q 

J = 0, Co = 


• + 


while max '^ub fo.C'j) ^ 




4 : 

5 : 

6 

7 


I^J+1 


|Aub(m;C'j)| 

Aub (Af;C'j)—A lb (Af;C'j) 
|Aub(a‘;C'j)| 


do 


arg max 
/iSS 

C'j+i ^ Cj u fij+i 

recompute Aub(i^;Cj+i) and Alb(/4;Cj+i) according to Q 
J ^ J+ 1 

end while 


2.3 Computational complexity 

Let us briefly summarize the computations performed by SCM. The bounding box B for 
y needs to be determined initially by computing the smallest and the largest eigenvalues 
of Ai, ...,Aq. Since each iteration requires the computation of the smallest eigenpair 
{Xi,Vi) of A{ij.i), this amounts to solving 2Q + J eigenproblems of size A x A in total. 
Verifying the accuracy of the current approximation on H and selecting the next param¬ 
eter sample requires computing Aub( 4 *; Cj) and Alb(i^; Cj) for all /r G H. In total, this 
amounts to solving J|H| LP problems with Q variables and at most 2Q + J constraints. 

2.4 Interpolation results 

As also discussed in [12], it is immediate to see that the bounds produced by SCM 
coincide with Amin(A(/x)) for all /r G Cj. The following theorem shows that the upper 
bounds also interpolate the derivatives of Amin(A(;u)) on Cj. 
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Theorem 2.1. Let Cj C D be finite and consider the upper bound Aub(/u;C'j) defined 
in Q. Given pLi G Cj in the interior of D, assume that 6i,...,6q : D —)• M are 
differentiable at pi and that Xi = Amin(^(^i)) is a simple eigenvalue of A{pi). Then 

VAuB ) C'j) = VAmin) 


with the gradient V with respect to p. 

Proof. Let Vi be an eigenvector associated with A* snch that ||nj ||2 = 1 and set yi := 
R{vi) G Tub(C' j)- By the definition Q, the relation 

Aub(m;C'j)= min e{p)'^ y = Oipf'yi (9) 

ye^uB (Cj) 


holds for p = Pi. The simplicity of Xi implies that yi is the unique minimizer. Combined 
with the facts that TubC^j) is a discrete set, pi is an interior point, and 6 {p) is continuous 
at Pi, this implies that ([^ also holds for all /r in a neighbourhood Ll C D around pi. 
Consequently, 


9Aub 

dp(Pi 


{tL Cj) 


de 

dp^Pl 


{Ti)'^yi, 


where p^Pi denotes the pth entry of p for p = 1,... ,P. 

On the other hand, the well-known expression for the derivative of a simple eigen¬ 
value [To] gives 


5A, 


dA 


dp^p) 


{A{pi)) = v*^^{pi)vi = v*[Y,j^){Ti)A, 


Q 


dOr. 


Q 


V'' t \ * /I ! \T 

= Vi, 


dp^P) 

q=l ^ 

5/iT) ' 


q I 


which completes the proof. 


□ 


As the following example shows, the result of 
lower bounds produced by SCM. 

Example 2.2. For p € D := [0,7r], let 


Theorem 2.1 does 


not hold for the 


A{p) = cos{p)Ai + sin(^)A 2 


cos{p) 


1 

0 


0 

-1 


-bsin(^) 


-1 

0 


It can be shown that y, the joint numerical range of Ai and A 2 , equals the unit circle 
around 0. Consider the sample set C 3 = {pi, p 2 , ps} = {0, |,7r}, with X mw iAjp^ )) = 
Amin(A(/X 2 )) = Amin(^(/.i 3 )) = —1- The resulting lower bound set TlbCC's) is the half¬ 
infinite box shown in Figure^ When minimizing 6{p)'^y for y G Tlb(C' 3 ), the minimum 
is attained at the vertex (-1,-1) for p G [7r/4,7r/2] and at the vertex (1,-1) for p G 
[7r/2, 37r/4]. Hence, 


Alb (it) 


— cos li — sin/X, /or G [7r/4,7r/2], 

cos IX —sin/i, /or/X G [7r/2, 37r/4], 
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Figure 2: Joint numerical range y (red circle) and lower bound set (yellow 

area) for the setting described in Example 


2.2 


yielding the following one-sided derivatives at /i = 7r/2.- 

'^LB((7r/2)“) = sin(7r/2) - cos(7r/2) = 1, 

^LB((7r/2) + ) = - sin(7r/2) - cos(7r/2)) = -1. 


In contrast, the exact eigenvalue is differentiable at 7r/2. Moreover, /2)) = 0 

is different from both one-sided derivatives o/Alb(^)- 


Theorem |2 .1 1 and Example 2.2 indicate that the lower bounds produced by SCM are 
asymptotically less accurate than the upper bounds. Indeed, this has been observed 
numerically m, implying a need to find more accurate lower bounds. However, the 
following theorem indicates that such an improvement is not possible without taking 
additional information on into account. 


Theorem 2.3. Let Cj = {/ii,... ,/ij} T D and consider the lower bounds XhBih] Cj) 
defined in 0 for a Hermitian matrix function A{/j) in affine linear decomposition 0. 
Let pL ^ D. If J < N then there exist matrices Ai ,..., Aq G , defining A{ix) = 

6i[pi)Ai + • • • + 9q{pl)Aq, such that 


Amin(^(/^)) — Alb(/^)C*j) and Aniin(^(h'*)) — Amin(^(/^i)) (fd) 


hold for i = 1,..., J. 

Proof. Let the columns of F G and V± G form orthonormal bases 

of V = spanjui,..., uj} and V'*', respectively, where each Vi denotes an eigenvector 
associated with Amin(^(/^i))- Moreover, let y-p G Tlb(C'j) C denote a minimizer 
of 0 for /i, that is, Alb(/^;C'j) = The rest of the proof consists of showing 

that the matrices defined by 


:= FFM,FF* + yji^qVxVf, qG{l,...,Q}, 


satisfy (10). 
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Given u G C^, we can write u = u\> + u± with ny G V and u_l G V"*". For any /r G -D, 
we therefore have 

_ ^ 

u*A{fj.)u = + yfi,q\\u±\\f) = ul,A{n)uv + 0{^i)'^yj: \\u_i\\l. (11) 

9=1 


For y = yi, this yields 

U A(^Hi)u > Aniin(^(/li)) ||llV II 2 ■^min(^(Al*)) ll^-L II 2 ~ '^min(^(/l*)) I|ll|l 2 ) 

where we used that G 3^lb(C'j) implies > Xmm{A{iJ.i)) . Since equality is 

attained for u = Vi, this implies the second equality in (10). 

To show the first equality, we first notice that the definition of Alb(/^; Cj) implies 

R{'uv)\Wv\\2 > Alb(/i;C' j)||'wv||2- 


Inserted into (11) for y = y, this yields 

u*A{y)u > Alb(ii;C' j)||iiv ||2 + Alb(ai;C'j)||ii ±||2 = Alb(/i;C'j)||ii|| 2 - 


Since equality is attained by any u G V"*-, 
completes the proof. 


this shows the first equality in (10) and thus 

□ 


Since the definition of the lower bounds in (fTl) only depends on 6 (u) and the eigen- 


are identical with those for A{iJ,). For A{y), the lower bound Alb(/i;C'j) coincides with 
the exact eigenvalue at an arbitrary fixed J1 ^ D. Hence, additional knowledge, beyond 
the eigenvalues at Hi, needs to be incorporated to improve the lower bounds. 


values at fj,i, the lower bounds for the matrix function ^(/r) constructed in Theorem 


2.3 


3 Subspace approach 


In this section, our new subspace approach is presented that takes eigenvector informa¬ 
tion across different parameter samples into account and offers the flexibility to incor¬ 
porate eigenvectors for larger eigenvalues as well. 

Given Cj = {yi,... ,/ij} C D, suppose that for each sample we have computed 
the i'>\ smallest eigenvalues 

Ai = aS^) < Af) < • • • < \f 

of A{yLi) along with an orthonormal basis of associated eigenvectors ..., G 

C^. To simplify notation, we assume that an equal number of eigenpairs has been 
computed for each /xi, ..., yj, although this is not necessary. The eigenvectors will be 
collected in the subspace 


V{Cj,t) := spanjyj^^ ... 


.W 


^(1) 


( 12 ) 


In the subsequent two sections, we discuss how the information in V{Cj,i) can be used 
to compute tighter bounds for Amin (H(y)). 
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3.1 Subspace approach for upper bounds 

Given the subspace V(C'j, tj from (12), we define an upper bound set analogously to ([^: 

3^sub(Cj,^) := {R{v) ■ V e V(Cj,^)}. 

The corresponding upper bound for fj, £ D is defined as 

Asub(/^;C'j,£) := min 9{fj.)'^y. 

y&yviCj/) 

Clearly, we have Tub(C' j) ^ Tsub(C'j,-^) ^ y and thus 

Aub(/^;C'j) > Asub(/^; > Xmm{A{fj,)). 

To evaluate Xs\jb{i^', Cj, i) , we first compute an orthonormal basis V G t^^Nxje q£ 
V{Cj,i) and obtain 

Xs\JB{f^] Cj,£) = min = min 9{y)'^R(yw) 

v£V(Cj,e) 

Il™ll2 = l 

= min 9i{ij.)w*V*AiVw + --- + 9q{ij.)w*V*AqVw 
Il>"ll2 = l 

= X^i49i{y)V*A,V + ... 9Q{y)V*AQV) = X^^{V*A{y)V). (13) 

Thus, the computation of AsuB(/i, Cj, £) requires the solution of an eigenvalue problem 
of size J^ X Ji, with Ji usually much smaller than N. 


3.2 Subspace approach for lower bounds 


We will use a perturbation result to turn the upper bound (13) into a lower bound 


XsLB{lJ-]Cj,i) for n £ D. For this purpose, we consider for some small integer r < Ji 
the r smallest eigenvalues 

AsuB(/i; Cj, i) = < a® < • • • < A?) 

of V* A{iJ,)V, along with the corresponding eigenvectors wi,... ,Wr £ Let U £ 
be an orthonormal basis of the subspace U(fi) spanned by the Ritz vectors: 

:= spanlTrci,..., Vwr}- 

Moreover, let ?7 _l £ be an orthonormal basis of U^{yL) and denote the eigen¬ 

values of U]_A{y)U± by 

x(i) ^ x(2) ^ ^ AN-r) 

The transformed matrix 

[U,U^]*A{y)[U,UA_] = 


U*A{fi)U U*A{fi)U± 
UlA{y)U UlAly)U± 


clearly has the same eigenvalues as A(y), while the perturbed matrix 


U*A{y)U 0 
0 UlA{y)U^ 
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has the eigenvalues { Ay ^, ■ ■ •, Ay ^} U { , ■ ■ ■, A^ }. Applying a perturbation result 

by Li and Li [22] to this situation yields the error bound 

\Xram{A{n)) 

with the residual norm 


p := \\UlA{p)Uh = \\A{p)U - U{U*A{p)U)h 

and the absolute gap 5 := |Ay^ — A^^||. Rearranging terms thus gives the lower bound 

2p2 


/(A[^x) < Amin(A(/i)), with f{jf) := min {X^^’,v) - 


( 1 ) 


|Ay^ - r?| + \/|Ay^ -r/P + V 

(14) 


This lower bound is not practical so far, as it involves the quantity X^l, which would 
require the solution of a large eigenvalue problem of size (N — r) x (N — r). 

Lemma 3.1. The function / : M - 
increasing. 

Proof. See Section [Aj 


defined in (14) is continuous and monotonically 


□ 


Lemma 


3.1 


implies that /(??) remains a lower bound as long as rj < A^|. To 


marize, our subspace-accelerated lower bound is defined as 


Xslb{p; Cj,£) := min {X^\v) - 


2p2 


1 -^?^ -v\+ V |Ay ^ + V 


sum- 


(15) 


for a lower bound rj of A^|. 


3.2.1 Determining a lower bonnd for A 


( 1 ) 

‘W-L 


The lower bound for X^l = Ai^ 
adapting the ideas from Section |2.1 
for Xmm{A{p)) by solving the LP 


{U^A{p)U±) needed in (15) will be determined by 


Let us recall that SCM determines a lower bound 


Alb(/^; Cj) = 


mm 

y^yLsiCj) 


OihVv, 


(16) 


with 3 ^lb(C'j) •= {y E B : 6{fj,i)'^y > A^, i = 1,..., J} and the bounding box B defined 
in To simplify the discussion, we always assume in the following that Tlb(C'j) is 
a simple polytope with no degenerate facets. Then there exists an optimizer G 


of (16) such that there are Q, among 2Q + J, linearly independent active constraints 


In other words, satisfies a linear system 


= V’, (17) 

where 0 G jg invertible and each equation corresponds either to a constraint of 

the form 9{pi)^y^ = Aj or to a box constraint. In the following, we tacitly assume that 
at least one of the active constraints is a non-box constraint. 
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Establishing a lower bound for is equivalent to determining rj such that rj < 
u*j_A{fj,)u± holds for every u± G U^{^) with ||u_l ||2 = 1- The restriction of m_l to a 
lower-dimensional subspace can be used to tighten the non-box constraints in ( |16[ ). 

Lemma 3.2. With the notation introduced above, let A* = diag(A^^\ ..., A^^^) and Vi = 
,..., vf'^]. If N — r > r then 


u]_A{ij,i)u± > Ai -b Pi, 


where Pi is the smallest eigenvalue of the matrix 

[Ai - Vk) - V*UU*Vi{Ai - 


Proof. Using the spectral decomposition of A{fii), the result follows from 


min u*j_A{fii)u± > 

uj _ GV(-^ 

l^±ll2 = l 


min ulViAiV*u± + Af+^^ul(/ - VV*)u± (18) 

Il“±ll2 = l 

-b min u*j_Vi{Ai - u± 

ll“± Il2=l 

+ Xrnin{UlVi{Ai - xf^^hi)V*Ux) 

+ Xmin{V*UxUlV{Ai - xf+^he)) 

+ X^iniih - V*UU*V){Ai - xf+^he)) 

Xi + X^i^iAi - Xih) - V*UU*V{Ai - xf+^he)), 


where we used in the third equality that the negative eigenvalues of the matrix product 
U^Vii^Ai — lijVfUx do not change under a cyclic permutation of its factors. □ 


Using the values of Pi dehned in Lemma 3.2 


m 


we update the right-hand side pj G 

(0 ^ follows: If the kth equation corresponds to a non-box constraint 9{^i)^y = A*, 
we set 'tpk '■= 'Pk + Pi = Xi + Pi and, otherwise, V’fc := 'Pk- Since 0 is invertible, the 
solution of the resulting LP 


inf 6{fj,)'^y subject to Qy > ip 

y 

is trivially given by 

:= 0“^^. (19) 

This finally yields the desired lower bound 

:= 0(//)% < aJJI = X^u.{UlA{y)Ux)■ 

Remark 3.3. The choiee ofr, the dimension of the Ritz subspaee U{y), requires some 
consideration. For r = 0, U±{y,) = yields no improvement: Aslb(/u; Cj, £) = 
X]_,^{y; Cj, t). Intuitively, choosing r = 1 will be most effective when the second smallest 
eigenvalue of A{y) is well separated from the smallest eigenvalue. Otherwise, one may 
benefit from choosing slightly larger values ofr. In practice, we choose r by taking the 
maximal value of Aslb(m; Cj, i) for a few small values of r = 0,1,2 ,.... 
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3.3 Algorithm and computational complexity 

The implementation of the proposed lower and upper bounds requires some care in order 
to avoid unnecessary that involve quantities of size N, the original size of the problem. 

Computation of p. The quantity p = \\A{p,)U — UAu\\2 with Ajj = U*A{p)U = 
diag(Ay\...,Af^) can be computed by solving an r x r eigenvalue problem: 

= \raU{A{p)U{p)-U{p)A{p)r{A{p)U{p)-U{p)A{p))) 

= XraUU{p)*A{p)*A{p)U{p) - A{pf). 

Computation of and U* A{p)* A{p)U . By the affine linear decomposition ([^ , 

v*A{p)v = eiip)v*AiV + • • • + eQ{p)v*AQV. 

A standard technique in RBM, we compute and store the Ji x Ji matrices V*AqV, 
and update them as new columns are added to V. In turn, the computation of 
V*A{p)V, which is needed to evaluate the upper bound for every p G E, be¬ 
comes negligible as long as Ji <C N. Similarly, the evaluation of U*A{p)*A{p)U 
needed for p becomes negligible after the precomputation of V*A*AqiV for all 
q,q' = 


Choice of i. Clearly, a larger choice of i can be expected to lead to better bounds. On 
the other hand, a larger value of i increases the computational cost. Intuitively, 
choosing i larger than one appears to be most beneficial when the gap between the 
smallest and second smallest eigenvalues is small or even vanishes. One could, for 
example, choose i such that — A^^^ exceeds a certain threshold. However, in 

the absence of a priori information on eigenvalue gaps, it might be wisest to simply 
choose f = 1 for all pi. 

Algorithm summarizes our proposed procedure for computing subspace lower and 
upper bounds. Similarly as SCM, the algorithm requires the solution of 2Q + J eigen¬ 
value problems of size N x N for determining the bounding box in the beginning and 
the smallest i+1 eigenpairs in each iteration. Clearly, the latter part will become more 
expensive if £ > 1. However, we expect that this increase can be mitigated significantly 
in practice by the use of block algorithms. More specifically, when using a block eigen¬ 
value solver such as LOBPCG m and efficient implementations of block matrix-vector 
products with the matrix A (and its preconditioner), the computation of i smallest 
eigenvalues will not be much more expensive as long as i remains modest. 

Computing Xs\JB{p]Cj,i) and XsLB{p',Cj,i) for all p G E amounts to solving J|H| 
LP problems with Q variables and 2Q + J constraints, as well as J|H| eigenproblems 
of size Ji X Ji. As long as Ji <C iV, these parts will be negligible, and the cost of 
Algorithms and will be approximately equal. 

3.4 Interpolation results 

By definition, we already know that the bounds from our subspace approach are never 
worse than the bounds produced by SCM: 

X\,b{A]Cj) < \sBB{p]Cj,i) < Amin(A(/i)) < \s\]b{P] Cj, i) < \jjb{j\Cj), (20) 
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Algorithm 2 Subspace SCM 


Input: Training set H, affine linear decomposition such that A(/i) = 0i(//)Ai + • • • + 
9q{^)Aq is Hermitian for every fj, £ E. Relative error tolerance escM- 
Output: Set Cj C E with corresponding eigenvalues and eigenvector basis V of 
ViCjJ), such that ^sub(m;Cj,Q-^lb(m;Cj /) ^ ^ 


2 

3 

4 : 

5 : 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 


Asub(m;C'j,£) 

compute Amin(^q); Amax(^g) ioi q = I,..., Q, defining B according to Q 

J = 0, Co = 0 

while max -^sub A)--^lb C,Cj/) ^ do 




^J+i 


arg max 
/iSS 


-^suB (a* ;C'j/) 

-^sub(m;Cj/)—Aslb(m;Cj/) 


AsuB(Ai;C’j/) 


compute smallest eigenpairs (Aj^J.^^, ..., (Aj:[_]^, of A(/ij+i) 

Cj+i <— CjU /rj+i 

update V*AqV and V*A*Aq^V for all g, g' = 1,..., Q 

for jj, £ E do 

compute Asub(; 4;Cj+i,£) = XmmjV* _ 

compute p = VAmax(C(^)*A(/i)M(/r)C(^) - A(/r)2) 

compute Ufj, = argmiUj^g^j^g^,^^^^) 6{p)'^y and updated according to (19) 
p{p) £- 0{pfy 


db At) 


compute Aslb(ii; Cj+i, according to (15) 

end for 
J ^ J+ 1 
end while 


with equality at p = pi £ Cj. Together with Theorem 2.1, these inequalities imply that 
our upper bounds also interpolate the derivatives at pi. 

Corollary 3.4. For any i > 1 and any pi £ Cj that satisfies the assumptions of 
Theorem \2. 1\ it holds that 

VAsUB(/ 4i; Cj,f) = S/Xram{A{pi)) 


with XsvB{Pi',Cj,i) defined as in (13). 


Proof. By the assumptions, pi is an interior point of D and (20) implies that the in¬ 


equality Amin(d.(/r)) < Asub(/i; Cj, £) < Aub(/i;Cj) holds for all p va a neighbourhood 
of Pi- Combined with the result VAmin(d.(/ij)) = VX\j-B{pi-, Cj) of Theorem 2.1, this 
implies VXsvBiPi] CjA) = VXm\n{A{pi)). □ 

In contrast to SCM, it turns out that the subspace lower bounds also interpolate the 
derivative of Amin(^(/A)) at p £ Cj. To show this, we need the following lemma. 


Lemma 3.5. Let pi £ Cj satisfy the assumptions of Theorem \2. 1[ For any e > 0, there 
is a neighbourhood Ll C D around pi such that 


|di - A{)^(/r)| < e, 

Xf^ - pip) < e, 

hold for all p £ XI, where Ay^(-) and pi-) are defined as in Section 


( 21 ) 

(22) 


3.2 
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Proof. By construction, Ay^(//j) = Aj and thus the continuity of the smallest eigenvalue 
implies that © holds for all ^ in some neighbourhood hli around //j. It remains to 
prove (22). 


In the LP ([^ for determining Cj), which is trivially given by Aj, the con¬ 

straint 9{^iY'y = \i is active. Since we assumed that (Vlb(C' j) is a simple polytope with 
no degenerate facets, the continuity of 9{yi) implies that this constraint remains active 
in a neighbourhood 0 . 2 '- 9{yi)'^y^ = A* for all y G 122 ) where y^ is a minimizer of ([^ for 
determining XhBifJ-] Cj). 


By (18), the value of j3i defined in Lemma 3.2 satisfies 


Pi = min 

1 .AJ_ (fJ.) 


ulViAiV*ui_ + Af+^^ul(/ - ViV*)ux - Xi. 


The eigenvector belonging to the eigenvalue A* = Amin(^(A*i)) is contained in U 
for y = yi. Once again the simplicity of A* implies that the angle between and 
U becomes arbitrarily small as y approaches /ij. Therefore, for any e > 0, there is a 
neighbourhood O 3 such that 

Pi > min ulViAiV*ui_ + xf+^^ul{I - ViV*)u± - A, - ^ = Af^ - A^ - 


In summary, the vector y dehned in (19) satishes 


— Xi + Pi > Xi + x\ ^ — Aj — - 



(23) 


By the invertibility of 0, the vector y^ remains bounded in the vicinity of yi. To¬ 
gether with the continuity of 9{y), this implies that there is a neighbourhood 04 such 
that 

\{9{y)-9{yi)fy^\<^, V/r G O 4 . 


Combined with (23), this yields 


7]{y) = 9{yfy > Af^ - e 


which establishes (22). Setting 12 = n 122 O ^^3 n 124 completes the proof. 


□ 


The following theorem establishes the Hermite interpolation property of the subspace 
lower bounds. 


Theorem 3.6. Let yi G Cj satisfy the assumptions of Theorem 2.1 and, additionally, 
suppose that r < £ and > A^^\ Then 

'^XsBB{hi\Cj,l) = VAmin(^(A‘i))- 


Proof. By Lemma 3.5 and the simplicity of Amin(^()^i))) there is do > 0 such that 
vih) ^ Xy\y) + do for p sufficiently close to pi. Hence, the subspace lower bound (15) 
is given by 

2 p 2 


Aslb(i^;C' j,£) = Xy\p) - 


d -|- y^d^ -|- 4 / 9 ^ 


(24) 
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with p = \\U-^A{p)U \\2 and <5 = |Ay^(//) — p{p)\ > So- By Corollary 


3.4 


we have 


= VAsUB(/^i; = VAmin(^(/^i))' 


Since 5 is bounded from below, the result follows from (24) if the gradient of is zero. 
The assumptions AJ ' > and r < i imply that the invariant subspace belonging 
to the r smallest eigenvalues of A{pi) is simple and contained in V. Hence, standard 
perturbation results for invariant subspaces [36] yield p = 0{\\p — p,i\\ 2 ) for /r —)• //j and 
therefore V/?^ = 0, which completes the proof. □ 


3.4.1 A priori convergence in the one-parameter case 

In the following, we analyse the convergence of our subspace bounds for a special case: 
We assume that A{p) depends analytically on one parameter p G [—1,1] and, moreover, 
the eigenvalue Amin(^(A^)) is simple for all p G [—1,1]. 

Let Eji denote the open elliptic disc with foci ±1 and the sum of its half axes equal 
to R. Under the above assumptions, there is > 1 such that the (suitably normalized) 
eigenvector v{p) belonging to Amin(^(A^)) admits an analytic extension v : Eji^ —)> C^; 
see, e.g., [151132j . Note that v can be chosen to have norm 1 on [—1,1], see |32l Theorem 
XII.4], but this is not the case on Ef{^ in general. Let Cj = {pi,..., pj} contain the 
Chebyshev nodes pi = cos{‘^^'k) and set Vi := v{pi). The corresponding vector-valued 
interpolating polynomial is given by 


= h{p)vi 4-h ^j{p)vj 


(25) 


with the Lagrange polynomials ti,...,lj : [—1,1] —)• M. By extending a standard in¬ 
terpolation error result [25l Corollary 6.6A] to vector-valued functions (see, e.g., [TSl 
Lemma 2.2] for a similar extension), one can show that 


„ , , , ,,, {RERr^)M 


(26) 


holds for any 1 < R < Rq, with M = sup ||u( 2 ;)|| 2 . This result is utilized in the proof 

z^Epi 

of the following theorem, which shows exponential convergence of our subspace bounds. 

Theorem 3.7. Under the setting described above, the subspace lower and upper bounds 
for i = r = 1 satisfy 


^SVB{h',Cj,l) - XjniniMu)) < Cu R 

Amin(^(Ai)) “ AsLB(/a; Cj, 1) < ClR~^'^ 


(27) 

(28) 


for every p G [—1,1], with constants Cu,Cl independent of J and p. 

Proof. For (■ = 1, the subspace used in our bounds takes the form V = spanjui,..., vj{. 
The interpolating polynomial defined in (25) clearly satisfies pj{p) G V and hence (26) 
yields the following bound on the angle between V and v{p): 


min Iju - u (^)||2 < R 

v&V 


-J 


(29) 
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The inequality (27) now follows from existing approximation results for Ritz values; see, 
e.g., [291 Chapter 11], 

can 


To prove (|28|), we first note that the arguments from the proof of Theorem 3.6 

2p2 


be utilized to show that 

Cj, 1) = Asub(m; Cj, 1) - 


6 + + 4p2 


for sufficiently large J, where <5 > (5o > 0 for some 6 o not depending on /x or J. Since 
r = 1, the quantity p coincides with the residual of the smallest Ritz vector of A(p,) with 
respect to V. By (29) and |29l Theorem 11.7.1], we have p < R~^, which completes the 
proof. □ 


The maximal value of the exponent R in (27)-(28) depends on the gap between the 


smallest and the second smallest eigenvalue and the variation of A{p). This is discussed 
in more detail for a special case in [H Section 2.3.1]. 

3.5 Residual-based lower bounds 


As we will see in the numerical experiments (especially in Example 4.5), the subspace 


lower bounds can sometimes converge rather slowly in the initial phase of the algorithm, 
in contrast to the subspace upper bounds. This slow convergence can be viewed as a 
price that needs to be paid in order maintain the reliability of the lower bounds. In 
the following, we will propose an alternative that is heuristic (i.e., its reliability is not 
guaranteed) and is observed to converge faster in the initial phase. 

The alternative consists of simply subtracting the residual norm from the upper 
bound: 

AsubCi^; Cj,tj - \\A{p)u - Asub(7‘; Cj,t)u\\ 2 , (30) 

where u with jjujj 2 = 1 is a Ritz vector belonging to the smallest Ritz value Asub(/^; Cj, 1) 
of A(/x) with respect to V. A basic first-order perturbation result for Hermitian matri¬ 
ces [29] implies that (30) constitutes a lower bound for an eigenvalue of A{p), but not 


necessarily the smallest one. There is a risk, especially in the very beginning, that (30) 
is actually larger than the smallest eigenvalue, see Section]^ for examples. However, in 
all numerical experiments we have observed that a small number of iterations suffices 


until (30) becomes a lower bound for the smallest eigenvalue. 


Remark 3.8. When using the residual-based lower bound (30), it makes sense to also 
adjust the error measure ([^ that drives the sampling strategy to 

\\A{p)u — Asub(7‘; Cj,tju \\2 
“If 1Asub(i^;C'j,£)1 ’ 

and stop the iteration when this error estimate drops below escM- 


4 Applications and numerical experiments 

In this section, we report on the performance of our proposed approach for a number of 
large-scale examples. Algorithms and have been implemented in Matlab Version 
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7.14.0.739 (R2012a) and all experiments have been performed on an Intel Xeon CPU 
E31225 with 4 cores, 3.1 GHz, and 8 GB RAM. 

We compare Algorithm with Algorithm by computing the maximum relative 
error ratio Q . Additionally, we compare the convergence of the bounds from Sections 
and towards the exact smallest eigenvalues by measuring the absolute error 

max |bound(//) -A„,in(A(/i))|, (31) 

for the corresponding bound, both with respect to the number of iterations and with 
respect to the execution time (in seconds). 

When implementing and testing Algorithms and we have made the following 
choices. We set the relative tolerance to egcM = lO”"^, the maximum number of iter¬ 
ations to Jmax = 200 and the surrogate set H to be a random subset of D containing 
1000 elements. The smallest eigenpairs of A(//j) have been computed using the Matlab 
built-in function eigs, which is based on ARPAGK [20], with the tolerance set to 10“®. 
The Matlab built-in function linprog with the tolerance set to 10“® has been used 
for solving all linear programs. In all experiments, we have used Algorithm with the 
number of smallest eigenpairs included in V set to i = 1 , since this already provided 
significant improvements over Algorithm For choosing r from Section |3.2[ we have 
tested all values r = 0,1,..., Q, see Remark |3.3[ 

Remark 4.1. A slight modification of Algorithm^ can significantly reduce the time 
spent on solving linear programs. For fj. G E, suppose that yLBih-) £ Tlb(C'j) is a 
minimizer of Q on Tlb(C'j). Let (Aj+i,nj+i) be the smallest eigenpair o/A(/rj+i) 
computed in iteration J + 1. If 6 {plY' yi,-Q{n) > Aj+i, we have yLBilfi ^ 3^lb(C'j+i) C 
Tlb(C'j), making also a minimizer of Q on Tlb(C'j+i). In this ease, we have 

Alb(/^; Gj+i) = XhBih^', Cj) and there is no need to solve the linear program in (0. 

4.1 Random matrices 

We first consider an academic example, where a random dense Hermitian matrix Ai G 

C^xTV 

is perturbed, to a certain extent, by random Hermitian matrices A 2 , ■ ■ ■ ,Aq G 

ij^TVxTV. 

A(/x) = Ai + P.2A1 p^qAq, 

where jj. = {p, 2 , • •., /Uq) G D = [0, 

Example 4.2. We consider Q = A, N = 1000, 5 = 0.2, with Ai, A 2 , A 3 , A 4 having real 
random entries from the unit normal distribution. The performances of both algorithms 
is shown in Figure^ The convergence of Algorithm^^fiattens after around 25 iterations 
and does not reach the desired tolerance, while the convergence of Algorithm 0 is mueh 
faster and reaehes the desired tolerance within fl iterations. We have also considered the 
optimized version of Algorithm\^ as described in Remark \f.l\ To make the comparison 
fair, we have compared it to a variant of Algorithm 0 where the subspace bounds are 
recomputed only when Cj) is recomputed. The influence of these modifications 

on the performances of both algorithms can be seen in Figure 0 The optimized version 
of Algorithm 0 requires a slightly larger number of iterations to converge to prescribed 
accuracy, but it still outperforms even the optimized version of Algorithm^ Since Al¬ 
gorithm^ converges quickly, there is no need to even consider the residual-based lower 
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bounds from Section 3^. but we still include the results in Figures and\^for the sake 
of completeness. Here and in the following, the star denotes the iteration from which on 
the residual-based lower bound (30) actually constitutes a lower bound for the smallest 
eigenvalue. 





(a) Convergence of the maxi¬ 
mum relative error ratio ([^. 


(b) Convergence of the er- (c) Convergence of the er¬ 
ror (311 for the bounds w.r.t. ror (31) for the bounds w.r.t. 


iteration. 


time. 


Figure 3: Convergence plots for Algorithms and applied to Example 


4.2 





(a) Convergence of the maxi¬ 
mum relative error ratio ([^. 


(b) Convergence of the er- (c) Convergence of the er¬ 
ror (311 for the bounds w.r.t. ror (31) for the bounds w.r.t. 


iteration. 


time. 


Figure 4: Convergence plots for the optimized versions of Algorithms and applied to 
Example 


4.2 


4.2 Estimation of the coercivity constant 


A posteriori error estimation in model order reduction techniques for parametrized 
PDEs, such as reduced basis method, requires reliable estimates for the coercivity con¬ 
stant [3l| defined as 


«(Ai) = inf 

u&X 


a{u, u] p) 


|2 

\x 


(32) 


Here, is the coercive symmetric bilinear form in the weak formulation of the 

underlying PDE on the domain H and X is the accompanying function space with 
the norm || • \\x induced by the scalar product {u,v)x = -|- a{u,v;Jl)j with 

T > 0 and a fixed parameter value Ji chosen to be the center oi D. A finite element 
discretization of (32) leads to the minimization problem 


a^ip) 


. v'^ A{pL)v 

v'^Xv ’ 


(33) 
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where A{fi),X G are the matrices discretizing a(-,-,/x) and (•g)x, respectively. 

Minimizing (33) is clearly equivalent to compnting the smallest eigenvalne of the gener¬ 


alized eigenvalne problem 

A{fj,)v = XXv. 

By compnting a (sparse) Cholesky factorization X = LL^, we transform g: 

L~^ A{p)L~'^w = \w. 


Hence, the matrices Ai appearing in Assumption 1.1 need to be replaced by 

L~^AiL~'^, i = 1,..., Q. 

In the following, we consider three numerical examples of this type from the rbMIT 
toolbox m- We only include brief explanations of the examples; more details can be 
found in m and 


Example 4.3. This example concerns a linear elastieity model of a parametrized body 
(see Figure 5a). The parameter pLi determines the width of the hole in the body while 


the parameter fi 2 determines its Poisson’s ratio. A discretization of the underlying PDF 
leads to the matrix A{p,) = Q = 16, ^ = (/ii,/X 2 ) and functions 

that arise from the parametrization of the geometry. We choose N = 2183 and 
D = [-0.1,0.!] X [0.2,0.3]. As can be seen from Figure^ The results are similar to those 
presented in Example \4.S\ with Algorithm^ converging in 31 iteration and Algorithm^ 
not reaching the desired tolerance. 


Example 4.4. This example concerns a stationary heat equation on a parametrized do¬ 
main (see Figure\6^. The parameter pLi determines the coefficient in the Robin boundary 
conditions while the parameter fi 2 determines the length of the domain. A discretiza¬ 
tion of the underlying PDF leads to the matrix A{pl) = Q = 3? 

/i = (/ii,^ 2 ) and functions Oi{pL) arising from the parametrization of the geometry and 
boundary conditions. We choose N = 1311 and D = [0.02,0.5] x [2,8]. As can he seen 
from Figure\^ the results are similar to those observed in Examples 4^ and 4-3. 


Example 4.5. This example concerns a stationary heat equation on a square domain 
divided into blocks (see Figure 7a). In each of the subdomains, one of the parameters 
/ii,..., /ig determines a coefficient of the PDF 


div 


1 

-fJ-i 



0 on Hi, i = 1,... ,9. 


A discretization of the PDF leads to the matrix A{p,) = 9i{fi)Ai, where Q = 10, p. = 
{pi,..., pq) and functions 9i{p) arising from the parametrization of the PDF coefficients. 
We choose N = 1056 and D = [0.1,0.5]®. As can be seen in Figure^ the performance 
of both Algorithms'^ and\^ is not satisfactory, due to the slow convergence of the SCM 
and subspace lower bounds. Only the subspace upper hounds converges at a satisfactory 
rate. In this example, the residual-based lower bounds clearly show their advantage. They 
become reliable after only 31 iterations. 
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( 0 , 1 ) 


( 1 , 1 ) 


(0, /ii + i) 


( 0 ,W - 5 ) 


(i/ii + 5) 


(il^l - 5) 


( 0 ,- 1 ) ( 1 ,- 1 ) 



(a) Geometry of the underlying PDE. 



(b) Convergence of the maximum relative error 
ratio 



(c) Convergence of the 
bounds w.r.t. iteration. 


error (31) for the 


(d) Convergence of 
bounds w.r.t. time. 


the error (31) for the 


Figure 5: Convergence plots for Algorithms and applied to Example 


4.3 


5 Extension to computation of singular values 


In Section 4.2 we have seen that the computation of coercivity constants can be formu¬ 
lated in terms of 0- For non-elliptic parametrized PDE one may have to resort to the 
inf-sup constant m defined as 


a{fi) = inf sup 


a{u, V] fi) 


(34) 


“sx rex IpIIxII^'IIx ' 

where a(-, •,//) is the bilinear form in the weak formulation of the underlying PDE and 
X. A finite element discretization of (34) leads to the minimization problem 

A{fi)v 


inf sup = inf sup 

ueSX re^'^ X vAXux Xv ygRiv 


L ^A{fi)L 

\\xh\\yh 


(35) 


where, once again, A{/j) and X = are the discretizations of and (•,-)v, 

respectively. Minimizing (35) is equivalent to solving the singular value problem 

o-min(A"^A(/x)L“^), 
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(a) Geometry of the underlying PDE. 


(b) Convergence of the maximum relative error 
ratio 



iterations time 


(c) Convergence of the 
bounds w.r.t. iteration. 


error (31) for the 


(d) Convergence of 
bounds w.r.t. time. 


the error (31) for the 


Figure 6: Convergence plots for Algorithm and applied to Example 


4.4 


which, in turn, is equivalent to computing 


(36) 


since (Tmin(A) = ■\/A min (A^A). The expression (36) can be recast in terms of ([^, with 
terms, with the matrices Ai^j and functions for i,j = 1,... ,Q defined as 


-1 ^Tv-Ia^-T 


Aij = L-^Aj X 


The SCM algorithm has already been applied to (36) but only with limited success, since 
having terms in the affine decomposition of A{fj,) further increases the computational 
cost by making the solution of the LP problem ([^ significantly harder. The faster 
convergence of the subspace-accelerated approach to ( [36| ) mitigates this cost to a certain 
extent. 
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(a) Geometry of the underlying PDE. 




(b) Convergence of the maximum relative error 
ratio 



(c) Convergence of the 
bounds w.r.t. iteration. 


error (31) for the 


(d) Convergence of 
bounds w.r.t. time. 


the error (31) for the 


Figure 7: Convergence plots for Algorithm and applied to Example 


4.5 


6 Conclusions 


Solving a parametrized Hermitian eigenvalue problem can be computationally very hard 
and SCM is the most commonly used existing approach. We have proposed a new 
subspace-accelerated approach, given in Algorithm As can be seen in Section |3.4t it 
has better theoretical properties than SCM. As can be seen in Section]^ it also improves 
significantly on SCM in practice, for a number of examples discussed in the literature, 
while having only slightly larger computational cost per iteration. For problems with 
small gaps between the smallest eigenvalues, as in Example 4.5 the convergence of 


the subspace lower bounds may still not be satisfactory. For such cases, we propose a 
heuristic approach using residual-based lower bounds. The proposed approach can be 
extended to the solution of parametrized singular value problems. 
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A 


Proof of Lemma 


3.1 


As a composition of continuous functions, the function / is clearly continuous. To prove 
monotonicity we distinguish two cases. First, let rj > Ay\ Then 

fiv) = Ay ^ - 2pV (v - Ay ^ +, 

which clearly increases as rj increases. Now, let rj < Ay\ Then 
/(r?) =r]- 2p^l (^A^^ - r/ + \j[j] - Ay^)2 + 

and 

fkn) = 1 - - - . 

(Ay ^ -r] + y (Ay ^ - ??)2 + y - 7/)2 + V 
Showing f'{rj) > 0, and thus establishing monotonicity, is equivalent to 

(A?^ - vf + V + (Ay ^ - ?/)\/(Ay^ -r?)2 +V > 2p^ 

(Ay ^ -ri)\J (A^^ - r/)2 + V > q >-(A^^^ - r/)^ - 
which is trivially satisfied for Ay ^ > r]. This completes the proof. 
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